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Abstract 

In the last few years we have been developing a Monte Carlo simula- 
tion method to cope with systems of many electrons and ions in the Born- 
Oppenheimer (BO) approximation, the Coupled Electron-Ion Monte Carlo 
Method (CEIMC). Electronic properties in CEIMC are computed by Quan- 
tum Monte Carlo (QMC) rather than by Density Functional Theory (DFT) 
based techniques. CEIMC can, in principle, overcome some of the limita- 
tions of the present DFT based ab initio dynamical methods. Application of 
the new method to high pressure metallic hydrogen has recently appeared. In 
this paper we present a new sampling algorithm that we have developed in the 
framework of the Reptation Quantum Monte Carlo (RQMC) method chosen 
to sample the electronic degrees of freedom, thereby improving its efficiency. 
Moreover, we show here that, at least for the case of metallic hydrogen, varia- 
tional estimates of the electronic energies lead to an accurate sampling of the 
proton degrees of freedom. 
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1 Introduction 



Modern theoretical methods in condensed matter physics and chemistry rely heav- 
ily on numerical simulations. The problem of solving the Schroedinger equation 
for many-body systems is too difficult to be addressed directly, even within the 
simplification provided by the Born-Oppenheimer approximation. In the most pop- 
ular practical approaches (Hartree-Fock (HF) and the Density Functional Theory 
(DFT) based methods Q]) the original problem is replaced by the problem of solving 
the time independent Schroedinger equation for a single electron in the field of the 
nuclei and the mean field generated by the other electrons. DFT is, in principle, an 
exact theory but the energy functional must be treated approximately for practical 
purposes. In the simplest Local Density Approximation (LDA), this exact theory 
becomes a self-consistent mean field theory. Extensions of LDA, such as Gener- 
alized Gradient Approximation (GGA) provide more accurate results but remain 
essentially at the level of an effective mean field treatment. Despite the mean field 
character, DFT schemes have proved to provide quite accurate results for many 
different systems 

In 1985, Car and Parrinello introduced an efficient method to couple standard 
Molecular Dynamics for classical nuclei with the electronic structure calculation at 
the level of LDA done "on the fly" to extract the nuclear forces [2] • Because the 
method allowed study of the statistical mechanics of classical nuclei with many 
body electronic interactions, it opened the way for the use of simulation methods 
for reaHstic systems with an accuracy well beyond the limits of effective force flelds 
available. In the last twenty years, the number of appHcations of the Car-Parrinello 
ab-initio molecular dynamics has ranged from simple covalent bonded solids, to 
high pressure physics, material science and biological systems. There have also been 
extensions of the original algorithm to simulate systems at constant temperature and 
constant pressurep], flnite temperature effects for the electrons |4||, and quantum 
nuclei f5'|. 

Despite recent progress, DFT suffers from well-known limitations, for example, 
excited state properties such as optical gap and spectra are less reliable. DFT 
shows serious deficiencies in describing van der Waals interactions, non-equilibrium 
geometries such as reaction barriers, systems with transition metals and/or cluster 



isomers with competing bonding patternspllHI- As a consequence, current ab-initio 
predictions of metallization transition at high pressures, or even prediction of phase 
transitions are often only qualitative. Hydrogen is an extreme case[3IHlini but even 
in silicon the diamond//3-tin transition pressure and the melting temperature are 
seriously underestimated [Tfl] . 

Another route to the ground state properties of a system of many electrons 
in presence of nuclei is the Quantum Monte Carlo methodfTTl Ifi], In its simplest 
form, an analytic many electron wave function is chosen on the basis of the vari- 
ational principle (Variational Monte Carlo, VMC) and the quantum averages are 
obtained by a Metropolis Monte Carlo simulation of the electronic coordinates. A 
more accurate representation of the ground state wave function can be obtained by 
projecting the variational wave function with the operator exp{—f3eH} where H 
is the many-body hamiltonian, and Pe is the projection time. Provided that the 
variational wave function is not orthogonal to the ground state wave function, the 
projected function tends exponentially fast to the ground state wave function as 
Pe CO. Since matrix elements of the above projection operator at large values 
of l3e are unknown for non trivial systems, a Trotter breakup in many (P) small 
imaginary time intervals (te = Pe/P) must be employed. In the configuration rep- 
resentation, each projection corresponds to a 3A^-dimensional integral which can 
be performed by MetropoHs Monte Carlo method provided that the propagator in 
imaginary time can be chosen real and can be interpreted as a probability distribu- 
tion. This is the essence of the Diffusion Monte Carlo method (DMC) which is an 
"exact" method for systems of bosons or boltzmannons. This means that all system- 
atic errors in a simulation are under control in the sense that they can be reduced 
as much as desired. Since electron are fermions, the above scheme fails because the 
imaginary time propagator must be completely antisymmetric under exchange of 
two electrons and therefore cannot be chosen strictly non-negative everywhere in 
configurational space. This is the origin of the infamous "fermion sign problem". 
In order to avoid the sign problem the "fixed node approximation" has been pro- 
posed and used routinely to perform fermion simulations [Sl. The energy calculated 
with this approximation is variational with respect to the position of the nodal sur- 
faces of the trial wave function. Over the years, the level of accuracy of the fixed 



node approximation for simple homogeneous systems, such as "^He and the electron 
gas, has been systematically improved by introducing more sophisticated nodal sur- 
faces (backfiow orbitals)^! EI- In more complex, inhomogeneous situations such 
as atoms, molecules and extended systems of electrons and nuclei, progress have 
been somewhat slower. Nonetheless, in most cases, fixed-node QMC methods have 
proved to be more accurate than mean field methods (HF and DFT)|lSj. Computing 
ionic forces with QMC to replace the DFT forces in the ab-initio MD, is more diffi- 
cult and a general and efficient algorithm is still missing. Moreover, the computer 
time required for a QMC estimate of the electronic energy is, in general, more than 
for a corresponding DFT-LDA calculation. These problems have seriously limited 
the development of an ab-initio simulation method based on the QMC solution of 
the electronic problem "on the fiy". 

In recent years, we have developed a different strategy based entirely on the 
Monte Carlo method both for solving the electronic problem and for sampling the 
ionic configuration SDace[THll5j. The new method, called the Coupled Electron- Ion 
Monte Carlo method (CEIMC) has been applied so far to high pressure metallic 
hydrogen where it has found quite different effects of temperature than CPMD 
based on the LDA forces JBl- Our present interpretation of the disagreement is that 
LDA provides a Born-Oppenheimer surface quite smoother than the more accurate 
QMC one and this strongly affects the structure of the protonic system at T > 0. 

The paper is organized as follows. The following section[2lis devoted to an outHne 
of the CEIMC method. We will not go into all details since two long articles have 
appeared on general aspects and early implementations of the method^! Ej- One 
of the new aspects that we have recently implemented in CEIMC, not described in 
those references, is the Reptation Quantum Monte Carlo projection of the electronic 
variational wave functional]. So in the snbsection l2.1l we review the RQMC method 
and in the following snbsection l2.2l we focus on the sampHng algorithm, we introduce 
our new scheme to improve efficiency and reliability of RQMC and we provide an 
analytical proof. In SectiorElwe report numerical results on the convergence of the 
new scheme with the projection time and with the Trotter time step. Finally, in 
sectional we conclude. 



2 The Coupled Electron-Ion Monte Carlo method 

CEIMC method is based on the Born-Oppenheimer (BO) separation between the 
slow nuclei and the fast electrons. This is in contrast with other Quantum Monte 
Carlo methods, Diffusion Monte Carlo fDMClflTl |H] or finite temperature Path 
Integral Monte Carlo fPIMCl [T8lll9j methods where electrons and ions are treated 
on the same footing. As usual, the BO approximation allows to overcome the 
limitations of the other QMC methods, while introducing an often negligible error. 

In CEIMC, the configurational space of the proton degrees of freedom at inverse 
temperature /? = (fesT)"^ is sampled with a Metropolis algorithm in which the 
difference between the BO energy of a proton state S and of a trial state S" is 
computed by an electronic ground state QMC calculation. The QMC estimate of 
the energy difference A = [E{S') — E{S)] has statistical noise which would bias the 
standard Metropolis algorithm. Unbiased sampling of the proton configurations is 
achieved by the penalty methodf^ which replaces the energy difference A in the 
acceptance formula by A + (/3crA)^/2, where a'^ is the variance of the energy differ- 
ence. Since ti^ > 0, the noise always causes extra rejections but this compensates 
for "uphill" moves accepted because of a favorable energy fluctuation. 

Several methods for computing energy differences in QMC are availahle. fTHIT^] . 
A simple and efficient method is to sample the electronic degrees of freedom from a 
distribution function which is the sum of the electronic distribution functions for the 
S and S' states (e. g. the sum of the the square of the trial wave functions in VMC). 
Averages of operators involving electronic degrees of freedom and a single proton 
configuration, say S, (for instance total energy, variance, etc.) are then computed 
by correlated samnling [TTl ITU [15] . For the typical size of the proton moves (between 
O.OlA and O.SA for classical protons depending on density and temperature) and 
the typical system size (up to 54 protons) we have investigated, this method is much 
more efficient than performing two independent electronic calculations for the state 
S and S' . 

In the ground state QMC methods, an electronic trial wave function must be 
chosen according to the physics of the system being studied. For the metallic phase 
of hydrogen, we have recently developed analytic functions which include backfiow 
and three-body correlations [23 • These wave functions are particularly appropriate 



to CEIMC since they have accurate energies already at the variational level, they 
have no adjustable parameters requiring optimization, and their computational cost 
is much less than using orbitals expanded in a plane wave basis typically used in 
QMC calculations [22]. 

In metallic systems, finite size effects are large and must be suitably treated. The 
common procedure is to repeat the calculation for systems of increasing size and 
extrapolate to the thermodynamic limit but this is impractical within CEIMC, since 
it would have to be performed for any proposed protonic step before its acceptance. 
A much better strategy is to use Twist Averaged Boundary Condition (TABC) 
|2^i[ll5) which reduces the finite size error in the energy to the classical 1 /N behavior. 
It consists in averaging the energy over the phase that the many body wave function 
can pick if a single electron wraps around the super-cell. This is equivalent to 
Brillouin zone sampling in the single electron approximation. Within CEIMC it 
does not cause a large increase in required CPU time/step. 

Finally a recent improvement of the method is the introduction of quantum 
effects for the protons, quite important in high pressure hydrogen. This is done by 
developing the thermal density matrix of protonic degrees of freedom on the BO 
surface in Feynman Path Integrals JHl E|- A similar technique in the context of 
Car-Parrinello method has appeared^T. We are not going to discuss the last two 
aspects of the CEIMC. While TABC implementation in CEIMC has been described 
in ref.^j, our implementation of PIMC for proton degrees of freedom in CEIMC 
will be the subject of a future publication. 

2.1 Reptation Quantum Monte Carlo Method 

To go beyond VMC electronic energies, we implemented a Reptation Quantum 
Monte Carlo algorithm (RQMC)^^' rather than Diffusion Monte Carlo algorithm 
(DMC). The implementation of the energy difference method is more straightfor- 
ward in RQMC, nor are averages of observables which do not commute with the 
hamiltonian biased. 

In RQMC the ground state wave function is obtained by constructing an imag- 
inary time path integral for the electronic degrees of freedom. If |5'o) is the trial 
state, the trial state projected in a "time" /?e/2, \^ ~ e^'^=^/^|^'o), will converge 



to the ground state for large (3e- Let us define the "partition" function 
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The energy is then defined as 
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where the averages of the local energy, El{R) = 3?(^'q ^(i?) il^'o(-R)), are with 
respect to the path average. In practice, the energy is computed as the average 
of the local energy at the two ends of the path. Here SR indicates the real part in 
the case that the trial function or Hamiltonian is complex. The energy E{l3e) is an 
upper bound to the fixed node energy for each value of /3e, it converges to this at 
large /3e, and its /?e derivative is strictly negative. This latter quantity is, in fact, 
minus the variance of the total energy 



The variance tends to zero for large enough values of (3e providing a useful signal 
for the convergence of the energy to the ground state. This is the zero variance 
theorem in RQMC. On the other hand, the variance of the local energy computed 
at either end of the path is the mixed estimator of DMC for a^{(3e). In practical 
implementations, it is desirable to keep /3e as small as possible to maximize the 
efficiency of the energy difference method. 

To compute the needed density matrix elements, we divide the projection time 
/3e into P time slices Te = /3e/P and make a semi-classical approximation for 
exp(— Te-ff). Our notation for a single electronic configuration is i? = {ri, . . . ,rjv}, 
while for the entire path is s = {Rq, Ri, . . . , Rp}- The probability distribution for 
a path is 



where U{R) = ^[In'^oiR)] and Ls{R,R') is the symmetrized link action for our 
approximation of the short time propagator. We have used the importance sampling 
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Green's function of the DMC propagator 
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where the force is F(R) 



V[/(i?) and A = h? /2me, which provides the symmetrized 
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An alternative form for the link action could be obtained through the pair action 
developed in finite temperature Path Integral MC[l^. However, we have not im- 
plemented this form and do not have a comparison of its efficiency. 

In order to impose the fixed phase constraint on the projected wave function, we 
must add to the link action a term of the form Lf^{R, R') — Are Jq drj \\7 (f>{X {r]))f 
where (l>{X) is the phase of the trial wave function at electronic position X and the 
integral is taken over all paths X{r]) with boundary conditions ^(0) = R, X{1) — 
R' . We have taken an end-point approximation for this term except for real wave 
functions in which case fixed-node boundary conditions were used. 

Note that in the expressions above, the dependence on the nuclear degrees of 
freedom was not shown even though all quantities depend on them. The probability 
distribution of an electronic path will be n(s, S) where 5" indicates the position of 
all nuclei. Because we have an explicit distribution of the electronic paths, it is 
straightforward to apply the importance sampling scheme for the energy differences 
by sampling the probability distribution [n(s, S) + 11(5, S')] where S and S' are the 
current and the trial protonic state, respectively. Note also, for VMC n(s, S) oc 
4*0(5, S)^ becomes the square of modulus of the trial wave function (no projection 
and i?o — Rp)- 

2.2 The "bounce" algorithm 

In the original work on RQMC[13, the electronic path space was sampled by a 
reptation algorithm, an algorithm introduced to sample the configurational space 
of Hnear polymer chains. The slithering snake or reptation method seems to have 



originated by Kron^| and by Wall and Mandel[2^. Given a path configuration s, 
a move is done in two stages. First one of the two ends (either Rq or Rp) is sampled 
with probability 1/2 to be the growth end Rg. Then a new point near the growth 
end is sampled from a Gaussian distribution with center at Rg + 2\TeF{Rg). In 
order to keep the number of links on the path length constant, the old tail position 
is discarded in the trial move. The move is accepted or rejected with the Metropolis 
formula based on the probability of a reverse move. For use in the following, let us 
define the direction variable d as c? = +1 for a head move {Rg = Rp), and d = —1 
for a tail move {Rg — Rq). In standard reptation, the direction d is chosen randomly 
at each attempted step. 

In the standard reptation algorithm, the transition probability P{s s') is the 
product of an attempt probability Td{s s') and an acceptance probability ad{s 
s'). Note that the path distribution given in Eq.l^J, because of the symmetrized 
link action does not depend on the direction d in which it was constructed. In the 
Metropolis algorithm, the acceptance probability for the attempted move is 
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n{s)Td{s ^ s') 

which ensures that the transition probability Pd{s — » s') satisfies detailed balance 

U{s)Pd{s^s')=U{.s')P^d{s' ~^s) (8) 

The autocorrelation time of this algorithm in Monte Carlo steps, that is the 
number of MC steps between two uncorrelated configurations, scales as [{Pe/Te)'^ /A], 
where A is the acceptance rate, an unfavorable scaling for large /3e- Moreover the 
occasional appearance of persistent configurations bouncing back and forth without 
really sampling the configuration space has been previously observed^^. These are 
two very unfavorable features, particularly in the present context, where we need 
to perform many different electronic calculations (at least one per protonic move). 
There is a premium for a reliable, efficient and robust algorithm. 

We have found that a minimal modification of the reptation algorithm solves 
both of these problems. The idea is to chose randomly the growth direction at the 
beginning of the Markov chain, and reverse the direction upon rejection only, the 



"bounce" algorithm. As far as we are aware, "bounce" dynamics has not been previ- 
ously investigated for RQMC, though Wall and MandeljJH mentioned it without a 
detailed proof and subsequent polymer simulations did not use bounce, perhaps be- 
cause the acceptance ratio in the polymer systems is much smaller than in RQMC. 
There is a related algorithm for directed loop algorithm on the lattice and for sim- 
ulations of trapped diffusion |28l I27j . 

What follows is the proof that the bounce algorithm samples the correct prob- 
ability distribution n(s). The variable d is no longer randomly sampled, but, 
as before, the appropriate move is sampled from the same Gaussian distribution 
Td{s s') and accepted according to the Eq. Q. To be able to use the techniques 
of Markov chains, we need to enlarge the state space with the direction variable d. 
In the enlarged configuration space {s,d}, let us define the transition probability 
P(s, d s', d') of the Markov chain. The algorithm is a Markov process in the 
extended path space, and it is ergodic as DMC method, hence, it must converge to 
a unique stationary state, T(s,d) satisfying the eigenvalue equation: 



We show that our desired probability n(s) is solution of this equation. Within the 
imposed rule not all transitions are allowed, but P{s, d s' , d') ^ for d = o?' and 
s ^ s' (accepted move), or d' — ~d and s ~ s' (rejected move) only. Without loss 
of generality let us assume d' = +1 since we have symmetry between ±1. Eq. Q 
with T{s,d) replaced by n(s) is 



Because of detailed balance Eq.JHJ, we have n(s)P(s, 1 — > s', 1) = W[s')P{s\ — 1 
s, — 1), which when substituted in this equation gives 




(9) 



n(s')P(s', -1 ^ s\ i) + Y. n(5)P(s, 1 ^ s', 1) = n(s'). 



n(s') p(s',-i^s',i) + 5]p(s',-i^s,-i) =n(s'). 
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Note that we have completed the sum over s with the term s = s' because its 
probability vanishes. The term in the bracket exhausts all possibilities for a move 



from the state (s', —1), thus it adds to one. Hence n(s) is a solution of eq. © and 
by the theory of Markov chains, it is the probabiHty distribution of the stationary 
state. 

3 Results 

In order to check the vaHdity of our proof we first appHed the bounce algorithm to 
an analytically solvable model, namely a one dimensional harmonic oscillator and 
obtained the expected results. For a realistic test, we compare the standard and 
the bounce algorithms for a fixed pair of protonic configuration {S, S') generated 
during a VMC run of liquid hydrogen at r^. = 1.31 and T — 5000K. We have 
considered Np = ^ 16 protons and electrons using analytic wave functions with 
3-body and backfiow terms at the T point (periodic boundary conditions) [21]. In 
the test, we fixed the electronic imaginary time step to — 0.04 /i^^ and the 
projection time to Pe = 0.2 h^^ which corresponds to 4 links. The key quantity in 
CEIMC is the correlation time tc in electronic MC steps of the energy difference A — 
Ebo{S') — Ebo{S) which determines, for a fixed length of the electronic run and 
for a given proton displacement, the noise level. The shorter the correlation time tc 
the larger the number of independent determinations of the energy difference. This 
implies smaller noise level and a larger acceptance for protonic moves, i.e. a higher 
efficiency of the algorithm. In fig. we compare the histogram of the correlation 
time tc of the energy difference obtained with standard reptation and with the 
bounce algorithm over 400 blocks of 10^ electronic steps. In both calculations the 
electronic acceptance rate is 0.89 but the noise level is 0.28 with standard reptation 
and only 0.14 with the bounce algorithm, in agreement with the observed correlation 
times. Note that, not only the average, but also the width of the distribution is 
roughly twice as large with standard reptation than with bounce dynamics. 

Next we study the convergence of the bounce algorithm with respect to Tg and /?£. 
We first consider protons on a bcc lattice to study the convergence of total energy 
and variance and to compare with DMC. As above, we consider Nc = Np = 16 at 
rg — 1.31, with the boundary condition 9 — 27r(0.4, 0.5, 0.6). Data obtained with 
runs of 10® electronic steps, are shown in fig. [2 At fixed /3 = O.IGH^^ we observed 
a roughly linear convergence (from below) of the total energy with Te (not shown) . 
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Figure 1: Histogram of the correlation time tc of the energy difference. Comparison 
between the reptation and the bounce algorithm for a path with 4 links. 

The results in fig. |2lare for Tg = 0.04i/~^ which may underestimate the energy by 
Q.imH I atom. Because of the high quality of the trial function, the ground state is 
reached with a very small projection time. Already at (3e — 0.6 the energy saturates 
at the value obtained with DMC (essentially infinite projection time, it is shown as 
a horizontal line in the upper left panel) . The remarkable linear dependence of the 
energy versus the variance below — 0.005 (upper middle panel) can be used to 
reliably extrapolate the energy to the — > cxd limit. 

In order to study the convergence of the energy difference and to estimate the 
scaling of tc with /?e, we consider a pair of successive protonic configurations for 
the same system generated during a CEIMC run at T=5000K, i. e. in the liquid 
state. At fixed /?e = we study the convergence with Tg in the range 

O.OIH^^ < Te < O.OSH^^, and at fixed Tg = 0.02H~^, we study the convergence 
with l3e in the range O.OSi/"^ < P < 9.6i7~^ which corresponds to 4 < P < 480 
time slices. In figure we show the first two moments tc and (7^ of high quality 
Gaussian fits to the histograms of tc- At fixed /3e, the rejection rate increases linearly 
with Te (not shown). However, successful moves are more effective and this results 
in the observed scaling ic ~ ~ t~°-^'^ (left panels) at least in the limited range 
of values of Te spanned. The behavior for increasing (3e at fixed Tc is more sluggish. 
Note that the rejection rate, 0.037 in the present case, does not depend on /3e. Both 
tc and o-g exhibit a somewhat erratic behavior but the overall scalings are quite 

2 

favorable. Note that Cc appears to scale roughly as tc . Although the quality of the 
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Figure 2: Bcc hydrogen at — 1.31. Convergence with (ie for the energy (left 
panel) and total energy versus the variance (right panel). In the left panel the 
curve is a shifted exponential fit while the horizontal dot-dashed lines represent the 
DMC result with its statistical error. . 

Gaussian fit remains good even at large /3e, for /3e > 1 the histogram of tc starts 
developing an small asymmetry with respect to the maximum with slower decay at 
large values of . 

In fig. 01 we report the related energy convergence study at fixed — 0.02iJ^^. 
In all panels, horizontal lines represent the variational estimate with its statistical 
error. In particular the panel a) shows that the energy difference AiJ/fesT used 
in CEIMC to perform the acceptance/rejection test is roughly independent of (5^ 
(neither is there any Te dependence at fixed /?e))- This result suggests that difference 
of the electronic energies at the variational level is accurate enough to perform 
CEIMC, at least in the present case of metaUic hydrogen with these analytical 
trial functions; we can sample the proton coordinates using VMC and compute 
the corrections to the energy and to the equation of state with RQMC for well 
equilibrated, statistically independent configurations. From panel c) we see that 
the projected energy is lower by b.TmH/ atom = 1%Q9K / atom with respect to the 
variational estimate, a significant change on the proton energy scale. In panels b) 
and c) are shown exponential fits to the data. Panel d) shows the energy vs the 
variance. As previously noticed, Hnear behavior is obtained for cr^ < 0.005. 

Finally, in order to test whether the VMC and RQMC computed BO surfaces 
have the same shape in the relevant part of the proton configurational space and 
not only at a single point, we have studied a system of Np = Ne ~ 54 atoms 
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Figure 3: Scaling of the average correlation time tc of the energy difference and 
of its variance cr^ for a fixed pair of proton configurations. Left panels show the 
behavior with Te at fixed /3e = 0.16H~^, while right panels show the /3e dependence 
at fixed Te = 0.02i/-i. 
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Figure 4: (3^ dependence of total energy, variance and energy difference for a pair 
of proton configurations {S,S'). The study is performed for = 0.02iJ^^. Dot- 
dashed lines represent the variational estimates with their error bars. In panel b) 
and c) the lines are exponential fits to data and in panel d) the continuous line is a 
linear fit in the region < 0.005. 



at = 1 and T = lOOO/C with zero twist phase (F point). Comparison of our 
VMC pair correlation functions with CPMD-LDA results[22| at this thermodynamic 
point has recently appeared 16j. The RQMC calculation has been performed with 
Te = 0.02_ff"i and (3e = l.OH^^ and provides an energy of -0.41114(8)if/atoTO to 
be compared with the variational estimate of — 0.4087(l)i?/atoTO. The computed 
pressure is 17.47(l)Mbars and VMC and RQMC estimates are in agreement within 
error bars. Average correlation time and variance of the energy difference are tc = 
7.1, al = 2.3 and = 16.5, cr^ = 26.5 for VMC and RQMC respectively. Therefore, 
going from VMC to RQMC with the same efficiency requires electronic runs between 
two and three times longer. 

4 Conclusions 

In conclusion, we have developed a new sampling algorithm for reptation Quantum 
Monte Carlo which we have shown to be more efficient than the standard sampling 
scheme and to have a favorable scaling with the projection (imaginary) time. This 
new scheme, which requires a minimal change of existing codes, allows one to sam- 
ple long electronic paths with a limited effort. We did not observe the occurrence 
of pathological situations previously reported with the standard scheme where the 
direction was resampled each move. We have implemented the new sampling al- 
gorithm in the CEIMC method and found that the correlation time of the energy 
difference for a given pair of protonic configurations grows like the projection time 
to the power 0.15. This means, in practice, that the noise level in CEIMC will get 
only moderately worse with increasing projection time, i.e. approaching the ground 
electronic state. More important, we have found that the difference in energy be- 
tween the two configurations is not sensitive to the projection time, suggesting that 
CEIMC sampling with VMC provides accurate dynamics. This conjecture has been 
verified for metallic hydrogen at a single thermodynamic point. 

An interesting question that remains unanswered is how general our conclusions 
are. Since the trial wave functions used in the present application are particularly 
accurate, which is not generally the case, caution must be exercised in applying the 
algorithms to cases where the accuracy of the trial function is unknown. 

Early aspects of the CEIMC algorithm were developed in collaboration with M. 
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